In a Remote Sensing Context

Maximum likelihood multivariate techniques are used to classify pixels from remotely sensed data into specific crop and other land cover types. This approach, although useful in estimating the acreage planted to specific crops, tends to create a visually speckled (or noise filled) image where misclassification and/or missing crop areas occur in a field. Past research into reducing this effect included two types of pixel smoothing algorithms, edge preserving and classified spatial, both of which use the spatial information contained in the data (Winings, Ozga, and Stakenborg (1990)). This paper considers two additional approaches which also exploit the spatial component of the pixel data in the classification process. Simple Augmentation (Switzer (1980)), and Contextual Classification (Owen (1984)) methods are explored. A comparison of all the techniques, applied to the 1996 Arkansas Landsat Thematic Mapper data, will be performed so as to determine their effect on the crop acreage estimation using regression estimation. INTRODUCTION The National Agricultural Statistics Service (NASS) of the United States Department of Agriculture produces state and national crop acreage sampled estimates through personal enumeration of farm information from randomly selected areas of land (called segments) in various land use strata. These strata are usually defined according to the percentage of cultivation in the given area. For example one strata might be >75% cultivated while another strata might be non-agricultural, or less than 15% cultivated. A direct expansion estimator based on enumerated segment data can be used to calculate crop acreage estimates. NASS uses remotely sensed data to improve these crop acreage estimates in selected states. Satellite imagery yields auxiliary information on a larger scale than the sample segments alone. This additional information can be used in a regression approach to create a more precise estimate of the crop acreage. Field boundaries determined in the sample segment enumeration process are training input into clustering algorithms for the remotely sensed data. Using the clustering results, we classify entire satellite scenes and then aggregate the data by land use strata. These aggregate values are placed into a regression estimator which yields our final estimate (Graham (1993)). Three measures of quality are used to determine the numerical accuracy of the remote sensing classification process. The standard remote sensing measures are percent correct and commission error; NASS also uses the regression coefficient of determination from the estimation process. Percent correct, for each ground cover, is the percentage of pixels from the fields of the segments correctly classified into that cover. Commission error, also derived from the ground truth, is a percentage of pixels classified as one cover but which actually belong to another cover, a type II error. The coefficient of determination, R, is a regression quality measurement. A regression is formed for each of the main crops, soybeans, rice and cotton. The explanatory variable is the number of classified pixels and the response is the crop acreage. A high coefficient of determination, its maximum being one, indicates that the number of pixels is accurately predicting the crop acreage. All processing for this evaluation was done using the PEDITOR system (Ozga and Craig (1995), Angelici, Slye, and Ozga (1986)) on the NASS MicroVAX and IBM-compatible PCs. New programs were added to PEDITOR to handle simple augmentation and contextual Bayesian classification methods. These programs run under Window NT which has fewer memory constraints for processing large arrays of numbers. The NASS acreage estimation process for Arkansas in 1996 created three analysis districts based on dates of satellite imagery. The test area chosen for this study was the Northern Analysis District from Landsat Satellite path 23 over Arkansas. This district is comprised of data from two satellite scene locations with two acquisition dates per location. Multi-temporal (in this case two date) data are used to improve the discrimination capabilities of the classification process: the first pass was acquired July 3, 996 and the second overhead pass on August 20, 1996. The digital imagery from each date was composed of three (of seven possible) bands from the Landsat Thematic Mapper: the visible red, near infra red, and mid-infra red. The tables of results in Appendix A show the percent correct, commission error, and R for the various crops and cover types in the district. Following Appendix A are five classified images of Craighead County, Arkansas (one county in the Northern Analysis District), plus five Zoomed-In images of a portion of eastern Jonesboro (the city in the center of Craighead county) and the neighboring cropland. There is one of each type of image for each of the methods explored in the paper. Note the legends on the county map apply to the Zoomed in image as well. THE STATISTICAL APPROACH Randomness is inherent in the remote sensing problem. The light reflected back in one part of a rice field will not correspond exactly to a light signature in another part of the same field. Hence for classification a statistical approach is required. The edge preserving smoother (EPS) and the simple augmentation (SA) are applied to the raw data, the raw data is then clustered, and classified using maximum likelihood estimation. The classified spatial smoothing (CSS) and contextual Bayesian classifier (CBC) require a classified image as input. In both cases the input was the Arkansas 1996 data classified by the maximum likelihood estimation methods. Since maximum likelihood estimation methods are used in all of the techniques, a brief explanation follows. Maximum likelihood estimation The data over a pixel centered at s = (u,v) is represented by the vector Z(s)  (Z1(s),  , Zp(s)), over some square lattice with sites s D where the Zi (s) is emitted or reflected energy in band i. The electromagnetic radiation is different for various surface types, and thus it is possible to classify surface covers. It is assumed there are {K: 1, ..., k} classes of covers, more specifically in our case, crops. Another assumption is that each pixel belongs to one and only one category. After assigning a cover to the pixels representing the surface, the area or acreage of each class is readily computed. The data Z(s) where s D, a square lattice relating to the entire satellite scene, belonging to the same class tend to form “clouds” of data points. The purpose of pattern classification is to partition the space  into decision regions, each region containing a specific class. A decision surface is the surface separating the regions. The pattern classifier then optimizes these decision surfaces to some criteria. In the maximum likelihood estimation classification procedure this criteria is to maximize the frequency of correct classification. (Swain (1973)). Let θ(s) represent the true, but unknown, ground class in which pixel s belongs. It is further assumed the data are independent and normally distributed by a conditional density f(|k) where μk are the means and the identical variance matrices are . Our criteria here is then to maximize the likelihood is μkZ(s) μkμk + logπ(k) where {π(k): k = 1, ..., K} are the prior probability distribution of θ(s). We determine μk and k by the sample-mean vectors and the sample-variance matrices computed from the training data. Recall that the purpose of the training data was that the true classes of the pixels were known (Cressie (1993)). The maximum likelihood estimation method performs well numerically. The percent correct for soybeans, cotton and rice are 72.35, 83.35, and 77.88 respectively. The percent commission error, a type II error, range from 13.71 for soybeans to 34.19 for cotton. The challenge for these alternative methods is either to statistically significantly increase the coefficient of determination for the three major crops, or to improve the visible image and not significantly decrease the R. The MLE has high R at 0.934, 0.97, and 0.968 for soybeans, cotton, and rice respectively. Visually, however, the MLE image is filled with too much speckling, or noise. Some of this noise may accurately be in the image but a large portion is due to the variability of the signatures. At a resolution of 30 meters, it is rare a mature crop field will contain several 30x30 meter bare soil, or non agriculture patches of land. The assumptions, that the pixels are iid and normally distributed, made by the MLE procedure are not realistic. “Spatial dependencies between pixels may be caused by scattering of reflected electromagnetic radiation from the surface of the earth, contamination resulting from over sampling and re-sampling, or spatial continuity of the ground classes” (Cressie (1993) pg. 503). In our case it is the later form of spatial dependency which we hope to take into account when classifying the data. REVIEW OF PAST TECHNIQUES Edge Preserving Edge preserving smoothing is applied to the raw data before the clustering or classification occurs. A moving five by five pixel window processes the center pixel. Within this five by five window, nine sub-windows each containing the center pixel are defined. They are similar to a “rotation bar around the center” (Winings, Ozga, and Stakenborg (1990), pg 648). The variances of the sub-windows are compared, the mean value of the sub-window with the least variance is the new value for the processed central pixel. This process is repeated for all pixels in a scene. “This technique assumes that variance is a measure of homogeneity” (Winings, Ozga, and Stakenborg (1990), pg 648). Hence the new value is relative to the sub-window containing the fewest edges which leads to the name. With multi-channel images the edge preserving algorithm uses the sub-window with the least total variance summed over all the channels (Stakenborg, (1985)). The algorithm is iterated by a user specified number of times. We iterated the procedure twice. Past experience showed two iterations to give the best results. (Winings, Ozga, and Stakenborg (1990). Small areas are removed during smoothing. Subwindows containing areas smaller than two pixels by two pixels of differing cover classes are considered to be mixed and not a good representation of the cover of interest. This filter’s unique properties are ideal because agricultural tracts are largely homogeneous. Stakenborg (1987) and Hill and Megier(1987) detail modified version of the edge preserving smoothing algorithm. The edge preserving smoothing increased the overall percent correct from 70.23 to 76.28. As to the individual covers, EPS improved in 13 out of the 14 categories. In 9 out of 14 categories the percent commission error was reduced. The R for soybeans, cotton and rice decreased slightly by -0.063, -0.008, and -0.043 respectively, non of which are statistically significant. The R for the maximum likelihood method are already high so there is little room for improvement. Visually, the classified image is noticeably more pleasing than the MLE image as a lot of the noise was eliminated. One thing to notice in the Zoomed-In image is that mean filters blur the edges between regions and are sensitive to outliers. Perhaps the median should be used. This would result in preserved edges which are not sensitive to outliers; however, median filters tend to round the corners of objects (Cressie (1993)). Most of the roads have disappeared in the EPS image, and there seems to be a lot of soybeans planted in the city--from ground truth we know this is not the case. The results from the edge preserving smoothing are listed in Appendix A under the heading EPS. Classified Spatial Smoothing Classified spatial smoothing is applied after MLE classification. Classified spatial smoothing is a method that utilizes a three by three pixels centered about pixel being classified. A dominant cover class for this window is determined by weights, and the center pixel is assigned to that cover. If the window is mostly one cover, areas of one or two pixels of a different cover are eliminated. “Areas of three pixels are reduced to one pixel. Effects on larger areas and areas with several covers are more complex.” (Winings, Ozga, and Stakenborg (1990)) The smoothing process always works to simplify the number of covers in the window. The overall percent correct went from 70.23 to 76.31. In 12 out of 14 crops the percent corrects rose. The percent commission error was reduced in every crop cover. The R increased slightly, but not significantly, in all cases: +0.005 for soybeans, +0.002 for cotton, and +0.005 for rice. Upon close inspection of the maps of Craighead county (the MLE map versus the CSS map), it can be seen that the CSS map contains less noise. From the Zoomed-In image we can see that the smaller roads have disappeared and there is no sign of field boundaries. Here however, the city is largely non-agriculture and woods, as it is supposed to be, as opposed to soybeans from the EPS image. REVIEW OF NEW TECHNIQUES Simple Augmentation One approach for the incorporation of spatial information into the classification of pixels is simple augmentation. Switzer(1980) considers a data vector Z(s), an N  1 vector, and augments it with the average of the four neighboring pixels to obtain a 2N  1 vector Z*(s0)  (Z(s), Z (s0)) where Z(s0) = (1/4){Z(u-1,v) + Z(u+1,v) + Z(u,v-1) + Z(u, v+1)}. This new augmented vector Z*(s0) is then classified using maximum likelihood estimation methods. This approach while easy to program can be difficult to implement due to the doubling of the size of the reflectance vector. PEDITOR is restricted by memory constraints to working with vectors no larger than 14  1. The overall percent correct rose from 70.23 to 73.34. This is not as high as either of the two previously discussed methods-but it is only slightly smaller. The percent correct for the crop covers rose in only 8 out of the 14 instances. SA’s most notable increase is that it correctly classifies 82.39 percent of soybeans. The R value for rice and cotton rose while soybeans declined; +0.008, +0.003, -0.019 respectively. Again none of these variations in the R are statistically significant. Looking at the Craighead county map it is readily apparent that there is less noise than the MLE map, yet the road structures and angularity of the individual fields are preserved. In the lower right corner of the blow-up you can even see boundaries around the rice fields in the middle of the cotton fields. We know these boundaries exist because they do not grow well together. The chemicals used on cotton kill the rice crop. Contextual Bayesian Classification In the contextual Bayesian classification method an effort is made to distinguish center pixels from boundary pixels. Switzer et al. (1981) found that for contiguous regions the information about the classification of the neighboring pixels, the four nearest neighbors, improved the classification of the center pixel. Let It denote the category of the central pixel t, and It*= (It , INt , IEt , ISt , Iwt) then the categories in the neighborhood. The categories’ model where in the plane is divided into convex regions by straight lines is considered. From this model P( It*| It) can be computed. A random process, the Poisson field, generates lines which partition the plane (Owen (1984)). Categories are assigned to each region, with neighboring regions being allowed to have the same classification. If the pattern of these regions is large in comparison to the size of the pixel then these neighborhoods can be described by one of three patterns (Cressie (1993)): X pattern T pattern L pattern i j j i i i i i i i i j i i i For simplicity each neighborhood is allowed a maximum of two categories, i and j. This model assumes the intensity parameter of the Poisson field, λ, is small compared to the size of the pixels. The following four local properties follow: (1) every neighborhood intersects one boundary line with probability β, and has negligible probability of intersecting two or more boundaries. The probability β was taken from a sampled portion of the classified image and was estimated to be 0.526 for the Northern Analysis Crop Reporting District. (2) Regardless of the intensity, the boundary in a given neighborhood is a line segmenthence only the Land Tpattern are possible. (3) For all neighborhoods, the conditional probability of an L-pattern given an intersection by a boundary is equal to α0.41. (4) Each rotation of the Tand Lpatterns are equally likely. Here, a more Bayesian approach is taken in that the joint prior probability distribution π(θ(u,v), θ(u1,v), θ(u,v1)) and the likelihood function f(z(u,v), z(u1,v), z(u,v1)|θ(u,v), θ(u1,v), θ(u,v1)) (conditional density also referred to as the likelihood function) are assumed known. Owen (1984) models the center pixel in reference to its nearest four-corner neighbors. His geometric probability model represents the distribution of classes in the plane. Under the assumption that θ is a stationary process the probabilities associated with these three patterns are as follows: P(X pattern) = 1β + π(i)β; P(T pattern) = π(j)(1α)β/4; P(L pattern) = π(j)αβ/4. From the given known distributions above the posterior distribution of the five levels can be computed from Baye’s Theorem, see Cressie (1993) for the formula. Assuming there are K classes the posterior distribution of θ(u,v) given the data {Z(u,v), Z(u1,v), Z(u,v1)} can be obtained, through summation, again see Cressie (1993) for details. The question now is how do the assumptions about the patterns fare in our real world example? In the northern portion of the analysis district 60.62% of the 36,943,642 pixels were found to belong to a pattern. The break up is as follows: 33.99% belonged to the X-pattern, 20.08% belonged to the T-pattern, and 6.54% belonged to the L-pattern. In the southern portion of the analysis district, out of a total of 8,468,739 pixels, 66.04% belonged to a pattern. The X-pattern contained 40.35%, the T-pattern 19.66%, and the Lpattern 6.03%. Although it would have been desirable for the overall percentage rate to be higher, our imposing the pattern structure on the data was not unreasonable. The CBC percent correct were higher in 5 out of 15 categories. The overall percent correct was 76.42 up 6.19 from the MLE classification. The percent commission error was lower in 10 out of 14 categories. The R was smaller for soybeans by 0.071, was larger for cotton by 0.012 and larger for rice by 0.004. Noting that none of these improvements are statistically significant, we cannot really determine a “best” spatial classification procedure. However, looking at the classified images one immediately notices the large amount of green (soybeans) represented in the image. It appears that if the classifier has a doubt it classifies the pixel as soybeans. This was due to the large prior probability soybeans had. A listing of the prior distribution is as follows: Corn 0.023; Waste 0.016; Win Wheat 0.01; Sorghum 0.01; Perm Past 0.02; Cotton 0.14; Soybeans 0.43; Woods 0.08; Water 0.01; Grass 0.00; Rice 0.2; Idle Crop 0.01; Non Agg 0.05; OtherC 10 0.00. If you look at the Zoomed-In image you can see the preponderance of green. The angularity of the fields is preserved, but any boundary around the fields is classified as soybeans. Again there are soybeans in the heart of the city of Jonesboro, which we know from ground truth to be false. This image is visually unsatisfactory. CONCLUSIONS Looking only at the coefficient of determination for the regression estimates, the winner for soybeans is the CSS method, for Cotton the CBC, and for rice the SA algorithm. The percent corrects show soybeans to favor CBC, rice EPS and cotton CSS. CSS held the lowest percent commission error for soybeans, and CBC held the lowest for both rice and cotton. What is interesting to note is that the Maximum Likelihood Estimation procedure didn’t win out once. Otherwise it is, at least numbers wise a tie. Visually however, it is quite a different story. Edge preserving eliminates all boundaries and roads and classifies soybeans in Jonesboro, AR. The CBC must also be thrown out for so heavily classifying soybeans into the picture. The MLE image is too noisy to be visually pleasing, although its numbers are okay. This leaves CSS and SA. If you look closely at CSS, you will see that the boundaries around the fields are either non-existent or are misclassified as some other crop rather than being classified as non agricultural. Very few roads are correctly classified. If you were going to overlay with a digital line graph the roads this would be a useful image. As a stand alone, visually, the SA is the better classification image. There is less noise than the MLE and the roads are correctly classified, as is the city. Field boundaries are in place and are not classified to other crop signatures. One could navigate by this image. Numerically, the difference between all the methods is not statistically significant. Visually, the SA stands out as the most pleasing and accurate. The one drawback to the SA method is that it doubles the parameter space. Currently PEDITOR can only handle 14 dimensions. A full satellite scene has seven channels. Thus it would only be possible to classify a uni-temporal full spectral scene, with PEDITOR as it exists now. This would be unsatisfactory. The second overpass is very important in maintaining the accuracy of our estimates with farmers intentions and early plant signatures being so similar. However, with the advent of Windows NT, it is hoped that as PEDITOR is moved and modified for this platform that the SA classification scheme can be implemented for full spectral multi-temporal data. APPENDIX A TABLE OF RESULTS